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1 Introduction 

This paper focuses on the accurate numerical rep- 
resentation of complex networks of evolving dis- 
continuities in solids, with particular emphasis on 
cracks. The limitation of the standard finite ele- 
ment method (FEM) in approximating discontin- 
uous solutions has motivated the development of 
re- meshing [1], smeared crack models [2, 3], the 
extended Finite Element Method (XFEM) [4-6] 
and the Phantom Node Method (PNM) [7]. 

We propose a new method which has some 
similarities to the PNM, but crucially: (i) does 
not introduce an error on the crack geometry 
when mapping to natural coordinates; (ii) does 
not require numerical integration over only part 
of a domain; (iii) can incorporate weak disconti- 
nuities and cohesive cracks more readily; (iv) is 
ideally suited for the representation of multiple 
and complex networks of (weak, strong and cohe- 
sive) discontinuities; (v) leads to the same solu- 
tion as a finite element mesh where the disconti- 
nuity is represented explicitly; and (vi) is concep- 
tually simpler than the PNM. 


2 Brief overview of the Phan- 
tom Node Method 

2.1 Introduction 

The static equilibrium of a body with volume U 
under body forces with density f (acting on fi) 
and traction t acting on the boundary Tq can be 
expressed in the weak form as: 

f e T (v) a (u) dU = [ v T f dfl + f v T t dr^ 

(i) 


where u is the displacement; v is the test function; 
e is the strain, related to u through the differ- 
ential operator relative to Cartesian coordinates 
4 as e = (u); and a is the stress, related 

to the strains through Hook’s law as cr = De, 
with D being the constitutive tensor. In the 
PNM, Fig. la, each real node i (characterised by 
its nodal coordinates and degrees of freedom, 
DoF, q^) is accompanied by a phantom node i' 
(with the same nodal coordinates x^, i.e. x^ = x*, 
and associated DoF q^/; in general, q^ ^ q i). 

2.2 Without a discontinuity 

If there is no discontinuity to be modelled by the 
element (e.g. before failure), the element is sim- 
ply a standard finite element. The vector of nodal 
coordinates is given as xq. In the case of Fig. la, 


Xj 7 = [xi,x 2 ,x 3 ,x 4 ] = [xi/,x 2 ',x 3 /,x 4 /] . (2) 


The Jacobian of the transformation from 
physical (x) to natural (£,) coordinates is: 


J 


dx _ dN 

dl ~ dT Xn ’ 


(3) 


where N is a standard matrix of shape functions. 
Assuming an isoparametric formulation, the dis- 
placement field u is related to the real DoF q 
through the same matrix N, i.e.: 


u = Nq. (4) 

In the case of Fig. la, q T = [qi, q 2 , q3, qj- 
The strains can also be expressed in terms of the 
real DoF q as: 

e = 4c(u) = ^(N)J- 1 q = Bq, (5) 
with B = 4 (N) J -1 , leading to: 

K = / B t DB det (J) dS, 


( 6 ) 


where F is the integration domain (in natural co- 
ordinates) corresponding to D (in physical coor- 
dinates), and to the vector of nodal forces: 

Q = / N T f det (J) dE + J N T t det (J) dr s 

(7) 

where r~ is the boundary corresponding to Tq. 
Equilibrium, Eq. 1, becomes then: 

Kq = Q. (8) 

Eq. 8 involves only the real DoF ; for complete- 
ness, the phantom DoF can be defined using con- 
straint equations so that they coincide with the 
real DoF; otherwise, they can be removed from 
the system of equations during assembly. 


2.3 With a discontinuity 

A discontinuity in the element can be predicted 
through generic stress-based criteria or otherwise; 
§3.7 delves in more detail into propagation mod- 
els and criteria. Once a discontinuity is predicted, 
the points at which it crosses the element are de- 
fined (points 5 and 6 in Fig. la). In this case, the 
vector of nodal coordinates is still given by 
Eq. 2, and the Jacobian J is still given by Eq. 3. 
However, the element is split in two partitions Da 
and Db, as indicated in Fig. la, and the displace- 
ment u is no longer related to the real DoF q by 
Eq. 4. Instead, the displacements ua and ub, in 
partitions Da and Db respectively, are interpo- 
lated separately from the respective DoF : 


u A = Nq A and u B = Nq B . (9) 


In the case of Fig. la, qj = [q B , q 2 ', q 3 , qj and 
qg = [qi,q 2 j q 3 ',q 4 ']- The strains then become: 

e A = (u A ) = jQ: (N) J -1 q A = Bq A (10) 
e B =£c(u b ) = J Q(N)J- 1 q B =Bq B .(ll) 


with B = /^(N)J 1 . Note that the B matrix 
(Eqs. 10 and 11) is the same for both partitions 
Da and Db- The stiffness matrices for partitions 
Da and Db are obtained by integrating over the 
corresponding domain in the natural space, i.e.: 


K a = f B t DB det (J) dE 
Je a 

K b = f B t DB det (J) dE. 


( 12 ) 


Note that the integrands for partitions Da and 
Db in Eq. 12 are the same; only the integration 
domains (Fa and Fb, respectively) are different. 
Similarly, the force vectors are: 

Q a = f N T f det (J) dE + 

+ [ N T t det (J) dr= and (13) 
J r = A 

Q b = f N T f det (J) dE + 

+ [ N T t det (J) dr=, (14) 

J ^B 

where again the integrands are the same for both 
partitions but the integration domains differ. Fi- 
nally, from Eq. 1, the equations of equilibrium are: 

K A q A = Q a and K B q B = Q B . (15) 


3 Floating Node Method 

3.1 Overview of the approach 

§ 2 shows that, in the PNM, phantom DoF do not 
need an associated position vector before the ele- 
ment is split (§ 2.2; in fact, before the element is 
split, the phantom nodes are also not needed, but 
their presence is convenient for implementation 
in existing FE codes). After the element is split, 
the nodal position associated with the phantom 
DoF is not the most suitable in terms of transfor- 
mation to the natural coordinate system nor in 
terms of integration. 

This observation forms the basis of the FNM, 
which we detail in this section. In the FNM, see 
Fig. lb, each real node i is characterised by its 
nodal coordinates x* and associated DoF q^; in 
addition, the element contains a suitable num- 
ber of floating DoF without pre-defined associated 
nodal position vectors. 


3.2 Without a discontinuity 

If there is no discontinuity to be modelled by the 
element (e.g. before failure), the formulation is 
the same as in the standard FEM or the PNM. 
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3.3 With a discontinuity 


Once a discontinuity in the element is predicted 
(see § 3.7), and thus the coordinates of the points 
which define the intersection between the discon- 
tinuity and the element defined, the element is 
split in two or more partitions. We will illustrate 
in detail firstly the most typical situation in which 
the element is split in two partitions, fl A and 
(Fig. lb). 

Unlike in the PNM, we define a vector of nodal 
coordinates for each partition fl A and Q B , x a and 
xb respectively. For the case in Fig. lb, these 
would be (c.f. Eq. 2): 

x n A = [ x 5j x 6 , x 3 , x 4 ] and x^ B = [x 4 , x 2 , x 6 , x 5 ] . 

(16) 

This is possible because X5 and X6 can be directly 
obtained once the failure criterion is fulfilled and 
the discontinuity is defined. Each partition has 
then a separate Jacobian (c.f. Eq. 3), 


Ja 


dx _ dN 
dt ~ d£ XfiA 


and 


dx dN 

Jb = « = « xn »' 
(17) 


The displacements ua and ub, in partitions 
fl A and n B respectively, are interpolated sep- 
arately from the respective DoF. Assuming an 
isoparametric formulation, 


u A = Nq A and u B = Nq B . (18) 


Note that the integrands for partitions fl A and Q B 
in Eq. 21 are different but the integration domain 
(£) is the same (the usual integration domain of 
a standard finite element). Similarly, the force 
vectors are (c.f. Eqs. 13 and 14): 

Q a = N T f det (J A ) dE + J N T t det (J A ) dr H 

Q B = ^N T f det ( J B ) dE + J N T t det ( J B ) dT s . 

(22) 

Finally, from Eq. 1, the equations of equilibrium 
are: 

K A qA = Qa and K B qB = Qb- (23) 

3.4 Different geometries for the discon- 
tinuities and integration 

Fig. lb (see also Fig. 2a) shows one particular 
case for the intersection between a quadrilateral 
element and a strong discontinuity. Another pos- 
sibility (needed also to represent accurately the 
crack tip, see § 3.6) is shown in Fig. 2b. Finally, 
other possibilities for the intersection may lead to 
different partitions of Q (Figs. 2c, 2d, 2e, 2f, 2g 
and 2h. Cohesive cracks (Fig. 2i) can be repre- 
sented by partitioning the integration domain Q 
into two quadrilaterals (fl A and Ub) and a surface 

OA). 


In the case of Fig. lb, which represents a strong 
discontinuity, q A = [q 5 , q 6 , q 3 , q4] and q B = 
[qi,q2,qe , ,q5 / ]- Note that there are four differ- 
ent sets of floating DoF: qs is different from q5/ 
and q6 is different from If, for instance, a 
weak discontinuity were to be modelled, only two 
sets of floating DoF would be included in the el- 
ement, and the DoF with a prime would coincide 
with those without a prime. 

The strains then become (c.f. Eqs. 10 and 11): 

e A = A (ua) = Q (N) J^qA = B A q A (19) 
£b = A (u B ) = ^ (N) J B x q B = B B q B (20) 

with B a = q (N) JX 1 and B b = q (N) Jg 1 . The 
stiffness matrices for partitions fl A and 0 B are 
(c.f. Eq. 12): 

K a = f BjDB A det(J A ) dE 

r ( 21 ) 

K B = I B b DB b det (J B ) dE. 


3.5 Representation of discontinuities 


Representing a discontinuity across an entire 
mesh implies passing information which defines 
this discontinuity to each element. With this in- 
formation, the element will follow the appropri- 
ate integration procedure (see Fig. lb). Level set 
methods are for instance well suited for this pur- 
pose. We present below an alternative to the level 
set method which offers some advantages for the 
representation of complex crack networks. 

We firstly define an edge status variable p, for 
each edge with global index j : 


m(j) 


— 1, if no discontinuity at edge E^; 

< 0, if E j is at a crack tip; 

1, if a crack crosses through E - . 

(24) 


and then define a dataset for each edge j : 


ma m ( 25 ) 
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where E,(j) contains the natural coordinates of 
the intersection of a discontinuity with edge j, 
if such intersection exists. The dataset is then 
made available to all elements that include edge 
j. The initial status of a discontinuity (prior to 
an analysis) can be defined through an array of 
initial values for the dataset in Eq. 25. 

3.6 Representation of a crack tip 

Considering Fig. 2j, suppose that the strong dis- 
continuity terminates at edge E 2 , i.e., the element 
to the right of E 2 has not failed. The existence of 
a crack tip at edge E 2 implies that q6 = qio- Clos- 
ing the crack tip can be achieved through inter- 
nal assembly of floating nodes or using constraint 
equations. 

Considering Figs. 3a and 3b, it is clear that 
closing the crack tip is not sufficient for an accu- 
rate representation of the crack tip; in fact; the 
transition from a strong discontinuity to no dis- 
continuity can lead to an artefact whereby the el- 
ement at the crack tip (see Figs. 3a and 3b) does 
not have the adequate topology for connecting to 
the adjacent nncracked element (lack of compat- 
ibility). Figs. 3c and 3d show that this can be 
mitigated to an extent using constraint equations 
which define the DoF at the crack tip through in- 
terpolation of suitable DoF of the element at the 
crack tip. 

A better solution is shown in Figs. 3e and 3f; 
the lack of compatibility can be avoided by us- 
ing a transition element at the crack tip (in fact, 
the code required for doing this is exactly the 
same that is used to represent weak discontinu- 
ities within the FNM). 

Fig. 3g shows that refinement elements can be 
used as well, leading to further improvements in 
the representation of the crack tip. If at least one 
refinement element is used (Fig. 3g), then the one 
step virtual crack closure technique (VCCT) can 
be used directly to model propagation of strong 
discontinuities, as detailed in § 3.7.2. 

3.7 Propagation of discontinuities 

3.7.1 Stress based criteria with cohesive 
formulations 

When a crack in a floating node element is rep- 
resented through a cohesive model, with the lat- 


ter having a stress-based failure initiation crite- 
rion embedded in its formulation, it is natural to 
also use an identical stress-based failure criterion 
to decide on the eventual propagation of a crack 
through an element. If the failure criterion for the 
cohesive element is: 

/(a) = 0, (26) 

then the propagation criterion (to activate the 
floating DoF so as to form a cohesive sub-element 
along the potential crack path) can be: 

/(a) = e, with — 1 < e < 0, (27) 

where 6 is a user-defined non-positive number 
such that the floating DoF are activated before 
propagation is due to occur. 

3.7.2 Virtual crack closure technique 

It is possible to propagate cracks without using 
stress measures directly; i.e. using fracture me- 
chanics concepts only. For this, Virtual Crack 
Closure Technique (VCCT) [8] is particularly well 
suited for using with the FNM. 

Consider Fig. 3g, where the element at the 
crack tip (refinement element labelled R) has not 
failed yet, but contains a weak discontinuity. The 
weak discontinuity is introduced to better repre- 
sent the crack tip. It is equivalent to a local mesh 
refinement, and can be extended to several ele- 
ments ahead of the crack tip. Let F be the inter- 
nal force vector at the crack tip and [q] be the 
displacement jump at the opposing edge in the 
wake element. Also, let Aw be the area of the 
crack surface in the element on the wake of the 
crack (for a 2 dimensional problem, Aw = Av^ 
where t?w is indicated in Fig. 3g and b is the thick- 
ness) and Act be the area of the crack surface in 
the element at the crack tip (for a 2 dimensional 
problem, Act = ^ct& 5 where £ct is indicated in 
Fig. 3g). Then, the energy release rates for modes 
I and II can be calculated as [9]: 

<28) 

<29) 

where F n and F t are, respectively, the components 
of F in the normal and tangential directions to 
the crack, and [g n ] and are respectively the 
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components of [q] in the normal and tangential 
directions to the crack. 

Once the energy release rates are known, a 
criterion of the form: 

f{Gi, Gu) = 0 (30) 

can be used to decide on the eventual propaga- 
tion of a crack through the element and on the 
direction of the crack. 

Calculating the energy release rates in Eqs. 28 
and 29 requires that the element at a crack tip has 
access to the displacement jump at the opposing 
edge in the element in the wake of the crack. This 
can be achieved by expanding the dataset for each 
edge j (Eq. 25) with more information to become: 

Mi)> Z(J), W(i)> n (i)> A (j) ( 31 ) 

where [qj (j'j is the displacement jump at the 
other (opposing) edge of the element in the wake 
of the crack (W in Fig. 3g), n (j) is the normal 

to the crack direction in element W, and is 

the crack area in element W. 

3.8 Implementation in existing FE 
codes 

The validation and application examples below 
were obtained through an implementation of the 
FNM in an User Element (UEL) subroutine in the 
commercial code Abaqus [10]. 

4 Validation 

4.1 Convergence and accuracy 

We compare the mesh convergence of the PNM 
and FNM in the evaluation of stress intensity fac- 
tors (SIF) for an edge crack propagating in mode I 
(Fig. 4a). The numerical evaluations for the FNM 
are performed with VCCT as presented in § 3.7.2, 
and for PNM with VCCT as implemented in the 
commercial software Abaqus [10]. The model has 
dimensions W = 9 mm, L = 37 mm and a = 
3 mm. The Young’s modulus is E = 200 GPa and 
the Poisson’s ratio is v = 0.3. The applied stress 
is a = 1 MPa. The plate is discretised uniformly 


with nw and til four-noded plane stress ele- 
ments, in the width and height directions, respec- 
tively. The four meshes created have = 

(9, 37); (18, 75); (36, 149) and (42, 173). The re- 
sults are summarized in Fig. 4b. The FNM can be 
seen to converge monotonically and more rapidly 
than the PNM used for comparison. 

We now evaluate the stress intensity factors 
(SIF) for a centre slant crack (Fig. 5a) obtained 
by the FNM against the corresponding analytical 
solutions [11] in mode I (K\) and mode II (An), 
for different orientations 0 of the crack. The nu- 
merical evaluations for the FNM are performed 
with VCCT as presented in § 3.7.2. For FNM 
with VCCT, when the crack separates the orig- 
inal element domain into a triangle and a pen- 
tagon, both the partitions shown in Fig. 2d and 
that shown in Fig. 2e are employed. The model 
has dimensions L = W = 10 mm and the hori- 
zontal projection of the crack is a cos 9 = 0.1 W. 
The Young’s modulus is E — 200 GPa and the 
Poisson’s ratio is v = 0.3. The applied stress is 
(7 = 1 MPa. The plate is discretised uniformly 
with a coarse mesh (80 elements across the width 
and 81 elements across the height). Plane stress 
conditions are assumed. The results are summa- 
rized in Fig. 5b. The data-points labelled Tnt. 1’ 
are obtained with FNM-VCCT using the parti- 
tion in Fig. 2d, and the data-points labelled Tnt. 
2’ are obtained with FNM-VCCT using the par- 
tition in Fig. 2e. 

4.2 Crack propagation 

A double cantelever beam (DCB) test is used to 
simulate a propagating crack for a case in which 
the analytical solution (using corrected beam the- 
ory [12]) is known. The material is representative 
of carbon/PEEK fibre reinforced composite [13] 
and is modelled as an orthotropic material. The 
geometry, boundary conditions and mesh are de- 
scribed in Fig. 6a; the width of the specimen is 
25.4 mm; the elements are 0.25 mm in length and 
the most refined elements near the mid-section are 
0.124 mm in height. The cohesive zone approach 
with a standard bi-linear law and a stress-based 
criterion as described in § 3.7.1 is employed to de- 
termine the initiation and propagation of a crack. 
Since the loading is purely in Mode I, the choice of 
failure initiation criterion and mixed-mode dam- 
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age propagation law are irrelevant; we therefore 
use the typical quadratic interaction criteria for 
both initiation and propagation. 

This case is analysed with both the FNM and 
the PNM (the latter implemented in the commer- 
cial software Abaqns [10]). For the FNM, the 
transition element shown in Fig. 3f is employed 
for the element in front of the crack tip. The par- 
titioning of the transition element into a floating 
node element containing a cohesive sub-element 
(Fig. 2i) is carried out as detailed in § 3.7.1 with 
e = 0. The propagation of the crack is then de- 
termined by the failure of the cohesive elements. 
The results are shown in Fig. 6b. 

5 Application: modelling of the 
growth of matrix crack den- 
sity in a cross-ply laminate 

In this section, we analyse the problem of inter- 
action between matrix cracks and delamination 
on a cross-ply [O 2 / 90d]s laminate of toughened 
glass/epoxy, tested in tension by Joffe and Varna 
[14]. In this problem, correctly capturing the ma- 
trix crack/delamination interaction mechanism is 
important for the accurate prediction of matrix 
crack saturation and consequent transition to de- 
lamination. 

Fang et al. [15] (see Fig. 7a) showed that us- 
ing non-matching meshes at a crack intersection 
(e.g. using the PNM for the 90° matrix cracks in 
the 90° ply and cohesive elements at the interface 
for the delamination) leads to an inaccurate repre- 
sentation of the displacement jump (and hence of 
the cohesive traction) at the interface. Capturing 
correctly the displacement jump requires further 
DoF at the intersection between cracks (Fig. 7b); 
the FNM method is particularly well suited to 
model intersecting cracks, capturing correctly the 
displacement jump at the interface. 

Based on the FNM, an element specifically de- 
signed for cross-ply laminates is formed with both 
real nodes and floating nodes, as shown in Fig. 7c. 
It makes use of the known position of the inter- 
face, so that the interface is not seeded with real 
nodes (Fig. 7d); instead, it is represented by cohe- 
sive elements formed with floating nodes. In this 
way, minimum seeding is required during prepro- 
cessing. 


Failure occurs essentially under tension (mode 
I) in the 90° ply and under shear (mode II) for 
the delamination. The choice of failure initia- 
tion criterion and mixed-mode damage propaga- 
tion law are therefore largely irrelevant in this 
case. We therefore use the typical quadratic in- 
teraction criteria for both initiation and propaga- 
tion. All material properties are taken from the 
literature [14, 16, 17]. Since the loading is uniform 
in tension, a 10% reduction on transverse tensile 
strength and mode I critical energy release rate 
is introduced in the element at the centre so as 
to initiate failure at the centre of the model. The 
model, Fig. 8a, represents half of the laminate, 
with symmetric boundary conditions applied on 
the bottom surface of the 90° plies. 

Figs. 8b and 8c show the failure pattern pre- 
dictions for the laminate when using the FNM el- 
ement from Figs. 7c and 7d. To demonstrate how 
crucial it is to have matching meshes at the crack 
intersections, a second model was created, differ- 
ing only in that only one cohesive sub-element 
is used to model the delamination in each FNM 
element (as in Fig. 7a), rather than two. This sec- 
ond model corresponds very closely to a model in 
which matrix cracks in the 90° ply are modelled 
with the PNM and the delamination is modelled 
independently using cohesive elements. The re- 
sulting crack pattern, at the same level of strain 
as in Figs. 8b and 8c, is shown in Figs. 8d and 8e. 
The second model predicts significantly less de- 
lamination than the FNM model. The simulation 
shows that delamination does not start from the 
element containing the matrix crack; instead, it 
occurs firstly in the elements next to the cracked 
element. This non-physical sequence of delamina- 
tion propagation is expected in the formulation of 
the non-matching meshes (Fig. 7a). 

Fig. 9 shows the crack density vs. applied 
strain predictions. While both models are able 
to capture the growth of crack density with ap- 
plied strain, only the first model (with matching 
mesh at the crack intersections) is able to predict 
saturation accurately; the second model contin- 
ues to predict an increase in crack density, albeit 
at a lower growth rate, after saturation should 
have occurred. This example thus demonstrates 
the capability of using the FNM to construct ele- 
ments for the modelling of specific geometries and 
complex crack networks. 
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6 Conclusion 

This paper proposes a floating node method 
which can be implemented in existing finite el- 
ement packages. The paper demonstrates that 
the floating node method has the following ad- 
vantages over alternative methods, in particular 
the phantom node method: (i) it does not intro- 
duce an error on the crack geometry when map- 
ping from physical to natural coordinates; (ii) the 
integration is simple, as it does not require numer- 
ical integration over only part of a domain; (iii) it 
leads to the same solution as a finite element mesh 
where the discontinuity is represented explicitly; 
(iv) it can incorporate weak discontinuities and 
cohesive cracks readily; (v) it can be readily com- 
bined with VCCT; (vi) it provides accurate pre- 
dictions for stress intensity factors under generic 
mode ratios; (vii) it is ideally suited for the rep- 
resentation of multiple and complex networks of 
(weak, strong and cohesive) discontinuities; and 
(viii) it can successfully predict certain interac- 
tions between matrix cracking and delamination. 
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(b) Floating Node Method 

Fig. 1. Comparison between the Phantom Node Method and the Floating Node Method. 


A FLOATING NODE METHOD 




(a) Strong discontinuity, 1 st case 
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(c) Strong discontinuity, 3 rd case 
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(e) Strong discontinuity, 5 th case 

A,_ 



(g) Intersecting cracks 



(i) Cohesive crack 



(b) Strong discontinuity, 2 nd case 



(d) Strong discontinuity, 4 th case 



(h) Weak discontinuity 



(j) Example of local DoF, vertices and edges numbering for 
the floating node element in Fig. 2g. To facilitate assembly, 
a pair of floating DoF are associated with each edge, and 
four floating DoF are internal. 


Fig. 2. Examples of different discontinuities that can be modelled by the Floating Node Method (see key 
in Fig. lb). 
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(a) Artefact in the transition from a 
strong discontinuity to no discontinu- 
ity. 


Element @ wake Element @ crack tip 




(b) Artefact when modelling the 
crack with cohesive elements. 


Element @ wake Transition element 



(d) The mid-point (or edge in 3D) of 
the cohesive element at the crack tip 
is constrained. 


(e) Using a transition element im- 
proves the representation of the crack 
tip. 


Element @ wake Element @ crack tip 


(W) (CT) 

AAA 


© e 

* © 


].y Interpolation 


; 

0 ( 

© 


(c) Artefact is mitigated by interpo- 
lating the crack tip’s DoF from neigh- 
bouring DoF. 
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(f) Transition element used with a co- 
hesive crack. 


Element @ wake Refinement element Transition element 


. (W) a (R) a (T) 

A A A A 



(g) Use of one refinement element improves representation of crack tip 
further, and allows direct use of VCCT. 

Fig. 3. Using local refinement elements and transition elements to represent the crack tip more accurately 
(see key in Fig. lb). 
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(a) Edge crack model 


(b) Mesh convergence of SIF 


Fig. 4. For this edge crack model, the FNM converges monotonically, unlike the PNM. 






A FLOATING NODE METHOD 


Fig. 5. 
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(a) Slanted crack model (b) SIF evaluations for different angles 6 

For this slant crack model, the FNM captures the SIF well in modes I and II for different angles 



(a) DCB model 
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(b) Predictions 


Fig. 6. DCB validation case, showing that, for the same mesh seeding, the FNM predicts accurately the 
force (P) vs. displacement curve while the PNM overpredicts the force. 



(a) Non-matching meshes at a crack intersection [15]. 
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(c) FNM element before matrix cracking. 



(b) FNM element representing intersecting cracks [15]. 
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(d) FNM element after matrix cracking. 


Fig. 7. Modelling the intersection between matrix cracks and delamination with non- matching meshes 
fails to capture the displacement jump; the FNM can address this. 
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(a) FNM model. Note that there is only one floating node element along the height (see Figs. 7c and 7d). 



(a) FNM model. Note that there is only one floating node element along the height (see Figs. 7c and 7d) 


(b) Accurate modelling of the transition from matrix cracking to delamination when using the FNM element from Figs. 7c 
and 7d. 



(c) Zoom of FNM mesh in Fig. 8b. 



(d) Model with non-matching mesh at the intersection, showing that some transition to delamination is not correctly 
captured. 



(e) Zoom of non-matching mesh in Fig. 8d. 

Fig. 8. Modelling the transition from matrix cracking to delamination in a cross-ply composite specimen. 
The red dots indicate failure at the corresponding integration points. 



Strain (%) 


Fig. 9. The saturation crack density is correctly captured using the FNM element from Figs. 7c and 7d. 
Non-matching mesh results in over-prediction of this density. Experimental data is from Joffe and Varna 
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